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Abstract 

In this paper we study the two dimensional motion of three linked rigid bodies 
moving through a fluid. The bodies change their orientation relative to each other 
in a way which mimics the swimming of fish. In contrast to previous simulations 
the bodies are connected by an elastic skin. The skin responds to the movement 
of the bodies and the pressure of the fluid and alters the flow around the bodies. 
In particular it prevents fluid moving between them. The system of bodies and 
skin is similar in appearance to a swimming leech or tadpole depending on the 
relative size of the bodies. We simulate the system using SPH, with three types of 
particles: liquid particles, boundary force particles determining the surface of the 
rigid bodies, and skin particles defining the elastic skin. The latter interact with 
each other, and with the boundary force particles to which they are anchored, by 
linear spring forces. The boundary force particles and the skin particles interact 
with the fluid particles by pair forces which are similar to the forces used in the 
Immersed Boundary method. The algorithm is based on a Lagrangian, and the 
equations of motion conserve linear momentum exactly and angular momentum 
very accurately. We compare the motion of rigid bodies with and without skin 
keeping the total mass of the bodies plus skin fixed. When the ellipses are identical, 
and the forward gait is used, the bodies swim ~ 13% faster when they are connected 
by skin, and they require less energy. When the ellipses have different sizes, with 
the front ellipse largest, they travel ~ 39% faster and use less energy. In the case 
of the turning gait, the identical ellipses turn faster with skin and use less energy, 
but the different sized ellipses turn more slowly with skin. The algorithm is simple 
and robust and can be applied to bodies of arbitrary shape and in domains which 
include free surfaces and stratified fluids. 



1 Introduction 



Most marine creatures swim by changes in body shape. During these changes the outer 
surface remains smooth because of the elastic properties of the skin and tissue of the body. 
In earlier papers we approximated the motion of marine creatures by considering linked 
rigid bodies moving in two dimensions in response to changes in the angles between them, 
but neglected the effects of skin and tissue (Kajtar and Monaghan 2008). The bodies 
we considered were three identical ellipses connected by virtual rods which allowed fluid 
to move between the bodies. Our high Reynolds number results were in good agreement 
with the two dimensional inviscid calculations of Kanso et al., (2005) and Melli et al., 
(2006) and for lower Reynolds numbers they were in good agreement with the viscous 
calculations of Eldredge (2006, 2007, 2008). 

Although this model of swimming creatures is very crude it gives a surprisingly 
accurate prediction of the motion of a leech (Kajtar and Monaghan, 2010). Nevertheless 
it is desirable to improve the model to bring it closer to the motion of actual marine 
creatures. In particular we wish to mimic the elastic properties of their bodies while 
eliminating the flow between them. A simple way to do this is to connect the bodies 
by an elastic skin. We do not claim that this is anything but a crude representation of 
actual tissue, but it represents important features of such tissue and opens the way to 
represent it more accurately. We note, in particular, that an elastic surface will deform 
under pressure forces from the liquid. 

We simulate the system using SPH, with separate particles for the liquid, the bound- 
ary of the rigid bodies, and the skin, and we derive the inviscid equations from a La- 
grangian variational principle. The viscous equations then follow by adding a standard 
SPH viscous term. We apply this algorithm to both straight line motion and to turning 
motion. The algorithm conserves linear momentum to within round-off error. The time 
stepping introduces relative errors in the conservation of angular momentum which are 
typically 10 -6 . However, because we approximate the infinite fluid by a periodic domain, 
there is a larger change in the angular momentum because periodic boundaries do not 
conserve the angular momentum of a particle system. An alternative approach to the 
simulation of swimming fish is the method described by Borazjani et al. (2008) and 
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Borazjani and Sotiropoulos (2008, 2009, 2010) who use an Immersed Boundary method 
(Peskin, 1977, 2002). In their method the fish body is triangulated and treated as an im- 
mersed boundary which moves in a specified way. This method is similar to our method 
except that we use particles for the entire system while they use particles only to specify 
the fish-body surface which is not allowed to deform under liquid pressure forces. Fur- 
thermore, their algorithm is currently designed for straight line motion so that the effect 
of torques on the body are not included. 

The plan of this paper is to discuss the SPH equations of motion and the modeling 
of the skin. We then compare the speed and power output of three identical linked ellipses 
with and without the skin both forward and for turning motions. Finally we apply the 
method to a system of three different ellipses linked as before. It is trivial to apply the 
algorithm to the swimming of bodies through free surfaces, and to studies of swimming 
in stratified media. 



2 SPH equations for the fluid 

The continuum equations we solve are the Navier-Stokes equations with boundaries formed 
by parts of rigid bodies and sections of skin. Apart from the introduction of the skin, the 
equations are the same as those we have simulated before (Kajtar and Monaghan, 2008). 
To simplify the paper we give the details of the SPH equations and refer the reader to 
the continuum equations described by Kajtar and Monaghan (2008). 

2.1 The acceleration equations 

In the following we use a and b for the labels of the liquid SPH particles, j for the label 
of boundary force particles on the rigid bodies, and a as the label for the skin particles. 
We write the equation of motion for the liquid particle a in the form 

^ = F a (flmd) + F a (body) + F a (skin), (2.1) 

where 

F a (fluid) = -J2 m b(^ + I ^ + Uab ) V " W " b > ( 2 ' 2 ) 

h \Pa Pb J 
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(p p \ Nb 
F a {body) = - m j( -j + ~k + J V ^ + m i r «i/(l r «il) 



(2.3) 



and 



F a (skin) = -J2 m i( I ^ + I -J + n -) VaW ™ + E E rn a r afT f(\r aa \). (2.4) 

\Pa Pa J k=laeS(k) 

F a (fluid) is the pressure and viscous force per unit mass due to the other fluid 
particles. F a (body) is the force per unit mass due to the rigid bodies. It consists of two 
parts. The first is a direct pressure interaction which is a result of deriving the equations 
of motion from a variational principle using the continuity equation as a constraint. The 
second is based on Sirovich's formulation of the effects of boundaries (Sirovich, 1967, 1968) 
which we have discussed elsewhere (Monaghan and Kajtar, 2009). Our prescription for 
this is similar to the Immersed Boundary method. A typical boundary particle j on the 
surface of the rigid body exerts a repulsive force m a rrij r j/(|r OJ -|) on fluid particle a along 
the line joining their centers. Here and elsewhere r a j = r a — Vj. Correspondingly, fluid 
particle a exerts an equal but opposite force mjm a Yj a f(\r a j\) The form of the function 
f(\r a j\) is chosen so that it mimics a delta function and provides a force on the fluid particle 
which is normal to the surface of the body to a very close approximation (Monaghan and 
Kajtar 2009). The force per unit mass due to the skin particles F a (skin) is identical 
except the summations are over skin particles. 

In these equations is the mass of particle b, Pb and pb are the pressure and density 
at the position of particle b. We use the same equation of state to determine P in terms 



of p as that used by Kajtar and Monaghan (2008). Further details are given in ^2.5 This 
equation of state makes the fluid weakly compressible. Il a b specifies the viscous interaction 
between particles a and b. We use the same form of the viscous interaction as in Kajtar 
and Monaghan (2008). W a b denotes the smoothing kernel W{r a — r 6 , h ab ) and V a denotes 
the gradient taken with respect to the coordinates of particle a. In this paper W is the 
fourth degree Wendland kernel for two dimensions (Wendland, 1995), and has support 
2h ah . In the present calculations the h a b used in W a b is an average h a b = (h a + h b )/2. The 
choice of h is discussed in detail by Monaghan (1992, 2005). In this paper we choose the 



initial h for any particle to be 1.5 times the initial particle spacing but, thereafter, it is 
determined by the local density. The total number of bodies is iV& and the total number 
of skin segments is N s . B(k) denotes the set of labels associated with body labelled k and 
S(k) denotes the set of labels associated with skin segment k. 

The acceleration of the center of mass of body k with mass takes the form 

where is a constraint force associated with the specification of the angles <p between 
the bodies. The torque equation is 

Ih ~w = E ( Tj ~ Rfc ) x fj + Tki 

where is the torque associated with the constraints. The constraint forces and torques 
are discussed further in §2.3[ 

The force f, on boundary particle j of body k is given by 



fj = -rrij m a (-£ + ^ + U aj J VjW a j + ^2m j m a v ja f(\v ja \) + fj(skin), 

\Pa Pj J 



(2.7) 



where the first term is the force on the body particle j due to the pressure and viscous 
stress of the fluid, the second term is the reaction force arising from the forces on the fluid 



due to the second term in (2.3). The third term is the skin force which we discuss in ^2.2 



The acceleration of skin particle a is due to a pressure interaction with the fluid, 
a repulsive force interaction with the fluid (these are similar to those discussed for rigid 
body boundary particles), and a force per unit mass due to neighbouring skin and/or 
body particles, 
d 2 r / p p \ 

-jjf = ~^2 m a ( -f + -f + n aCT J V a W aa + ^2m a r aa f(\r aa \) + f a (body) + i a {skin). 

(2.8) 

The third term is the interaction with the bodies to which the skin is anchored (§ 2.2). 
The interaction between any pair of SPH particles is along the line of centres, and the 
force on one particle is opposite to the force on the other. As a consequence, linear and 
angular momentum are conserved. 
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The force function between the skin and either the boundary or skin particles has 
the following form: 



where B = V 2 /rh, m is the mass of a liquid particle, and V is the maximum speed of the 
fluid. V is estimated at the beginning of the calculation and thereafter held constant. W 
is the Wendland cubic kernel normalized to 1 at the center. With q = r/h it has the form 



In the present two dimensional study a section of the elastic skin is a single line of skin 
particles connected by spring-like forces. Only neighboring particles interact. When the 
skin is stretched to length L, the skin tension is T = kL, where k is a constant spring 
force per unit length. The skin thickness is £ s , and the mass of each skin particle o is 
tti = p a £ s A s , where p a denotes the skin density (which is constant), and A s the initial 
skin particle spacing. 

Skin particles interact with the fluid with the same boundary force as body particles. 
Skin particles can move in response to the forces acting on them, whereas body particles 
only move when the body to which they are attached moves. A typical configuration of 
body, skin and liquid particles is illustrated in Figure [TJ 

The elastic force on a skin particle with label o due to a neighbouring skin particle 
with label a', is given by 



The skin particles are always labelled such that a + 1, a and a — 1 are contiguous. In the 
case where the neighbouring particle is a body particle (to which the skin is attached), 
the label a' is replaced by that of the body particle. 

The continuum limit of our skin shows that the speed of a transverse wave propa- 




(2.9) 




(2.10) 



2.2 Skin 



f CT = k(iv - r CT ). 



(2.11) 
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Figure 1: An illustration of the placement of body and skin particles. Filled circles represent body 
particles and are placed around the perimeter of the elliptical bodies. Shaded circles represent skin 
particles. In this case, the skin is taut, but in general it flexes under the pressure of the fluid. The open 
circles represent fluid particles. The following vectors along which the forces are calculated have also 
been indicated: on fluid particle a due to body particle j, and on skin particle a due to j. 



gating along the skin is 



v s = \f^. (2.12) 
V m 

The parameter k — T/A s is then given by 



vim 



«=-fcr- (2-13) 

We choose the skin parameters so that v s is comparable to the speed of sound c s of our 
slightly compressible fluid to ensure the CFL condition from both speeds is similar. The 
details are discussed in connection with time stepping. 

The form of body particle-skin interaction fj(skin) in (2.7) is determined by the 
fact that only one of the boundary particles on a given rigid body can connect with 
a specified section of skin. The first and last particles on each section connect with 
a boundary particle of a body. In general, for this two dimensional problem, a body 
has four such connecting boundary particles, while the first and last bodies have two 
connecting particles. For any given boundary particle of a rigid body it either connects 
to a skin section or it doesn't. If it does connect it does so by an elastic force term. Thus, 
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for (2.7) 

{ft(r CT — Tj) ; j connected to a, 
(2.14) 
0, ; j not connected to a. 

2.3 The constraints 

The angle 9 k which fixes the rotation of body k is defined as the positive rotation of a 
line fixed in the body from the x axis of a cartesian coordinate system fixed in space. For 
simplicity we assume the line fixed in the body is an axis of symmetry. The constraint 
conditions on the angles are 

V?m = @m+l — 9 m , (2-15) 

where m is the link number and (p m is a specified function. The form of the (p m determines 
the gait of the bodies. For the examples we consider here there are three bodies and two 
links as shown in Figure [2] with the skin removed for clarity. In the simplest case ip m is a 
function of t but, in general, it depends on other variables. For example, in a biological 
problem, it could depend on the centre of mass coordinates in such a way that the fish 
slows down when it enters a region where food is abundant. 

In addition to the constraints on the angles there are constraints associated with the 
links. We assume the link, or pivot, is at a distance 4 from the centre of mass of body k. 
The condition on the X components of the centres of mass of bodies k and k + 1 is that 
the X coordinate of the link between them is given by 

X k - 4 cos (6 k ) - X k+1 - 4 +1 cos (0 fe+ i) = 0. (2.16) 

Similarly, the Y constraint is 

Y k - 4 sin (9 k ) - Y k+1 - 4+i sin (0 fc+1 ) = 0. (2.17) 

These constraints enable the coordinates of the centers of mass of the bodies, and their 
angles 9 to be written in terms of those of any selected body. Similarly, by differentiating 
the constraint conditions with respect to time, the velocities X and Y and angular velocity 
Q of the bodies can be written as functions of the same selected body. The number of 
degrees of freedom (coordinates and velocities) of iV linked bodies in two dimensions is 
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therefore 6 compared with the 6N degrees of freedom of N independent bodies in two 
dimensions. If the f m are functions of t alone it is possible to reduce the equations of 
motion to those involving the coordinates and velocities of one of the bodies. This can also 
be done when the <p m are functions of both coordinates and time but it is inconvenient to 
eliminate variables and, in our view, simpler to take account of the constraints by using 
Lagrange multipliers. For that reason we use Lagrange multipliers even though, in the 
applications to be described in this paper, the ip m are functions of t only. For the case of 
three bodies we have two links and therefore 6 constraints. 




Figure 2: The configuration of the bodies (assumed to be ellipses). The skin is not shown. The link 
position is denoted by a filled circle. The straight line through a body passes through its centre of mass 
and is assumed to be rigidly attached to the body with one end attached to the link. The angles 9 are 
defined relative to a fixed direction in space (shown by parallel dotted lines) which is taken to be the 
x axis of a cartesian coordinate system in our calculations. The angles ip determine the gait and are 
specified functions of time. 



We denote the Lagrange multipliers for the X, Y and 6 constraints of link m by 
X^\ an d X^ respectively. Using standard methods for holonomic constraints (e.g. 
Landau and Lifshitz, 1976) we find the following expressions for the constraint forces F& 
and torques for the various bodies. For bodies 1, 2 and 3 respectively, 

F 1 = (\$\\$>), (2.18) 
F 2 = (-A?,-A?) + (Ag ) ,Ag ) ), (2.19) 
F 3 = (-A? ) ,-Ag ) ). (2.20) 

These constraint forces do not affect the total linear momentum of the bodies because 
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they sum to zero. 

The constraint torques on bodies 1, 2 and 3 respectively are 

Tl = -Aj 1} + X^h sin (0i) - Xy^h cos (flj), (2.21) 

r 2 = - \jp + X^l 2 sin (0 2 ) + A^ 2 sin (0 2 ) - (A^ 1} + X^)£ 2 cos (0 2 ), (2.22) 

r 3 = A^ + A^ 3 sin (0 3 ) - A^ 3 cos (0 3 ). (2.23) 

The Lagrange multipliers can be calculated quickly using a Newton- Raphson method. 
For example, with three bodies, the computational time is ~ 1CT 3 of the total compu- 
tational time. The details are given by Kajtar and Monaghan (2008). Extending the 
algorithm to 4 or more bodies is straightforward. 



2.4 The rate of change of density 

The rate of change of density of the fluid particles is 

d ^ = Y j m a w ail -V a W am (2.24) 

v 

where the summation is over the labels of all the fluid, rigid boundary and skin particles. 
In some formulations of SPH the summation is only over the fluid particles, but a better 
estimate of the velocity divergence, and therefore the rate of change of the density, is 
obtained by including the velocity of the boundary and skin particles. As mentioned 
earlier the inviscid fluid equations can be obtained from a variational principle using the 
continuity equation as a constraint. As a result the pressure terms in the acceleration 
equation of the fluid then involve all the particles. The density of the boundary force 
particles and the skin particles is kept fixed. In practice the changes in density are small, 
but they need to be correctly calculated to ensure that the pressure is estimated accurately. 

The h associated with any fluid particle can be obtained from h oc 1/p 1 ^ 2 though 
we calculate it in step with the density from 

^ = - A^. (2.25) 
dt 2pdt K ' 
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2.5 Equation of state and viscosity 

The fluid is assumed to be slightly compressible with an equation of state given by 

(2.26) 

where po is the reference density of the fluid. To ensure the flow has a sufficiently low 
Mach number to approximate a constant density fluid accurately, we determine the speed 
of sound by c a ~ 10V where V is the maximum speed of the fluid relative to the bodies. 
In this paper we take V = 2au where a is the semi major axis of the ellipse and uj is 
the frequency of the oscillation associated with the gait. The pressures of the body force 
particles and the skin particles are set to zero. 

The viscosity is determined by IT a b for which we choose the form (Monaghan 1997, 

2005) 

Pab\ r ab\ 

In this expression a is a constant, and the notation v afe = v a — v 6 is used. p ab denotes the 
average density \(p a + p b ). We take the signal velocity to be 

v si9 = \{c a + c b )(l + a) - 2^^, (2.28) 

where c a is the speed of sound at particle a (Monaghan 1997, although here we take v S i g 
to be half used in that paper and a is therefore a factor 2 larger). The kinematic viscosity 
can be estimated by taking the continuum limit which is equivalent to letting the number 
of particles go to infinity while keeping the resolution length h constant. By a calculation 
similar to that in Monaghan (2005) it is found for the Wendland kernel that the kinematic 
viscosity is 

v = ^Oihv sig . (2.29) 

o 

SPH calculations for shear flow agree very closely with theoretical results using this kine- 
matic viscosity (Monaghan 2006). Using these results we can write (2.27) in the following 
form 

Il ab = (2.30) 
If desired v can be replaced by using the Reynolds number. 
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2.6 Motion of the particles 

The position of any fluid or skin particle is found by integrating 

dr 

*= v - < 2 - 31 > 

The motion of a boundary particle can be determined from the motion of centre of 
mass and the rotation about the centre of mass. Thus for particle j on body k, 

dr- 

^=V, + Q,zx(r 3 -R t ), (2.32) 

where, in this two dimensional problem, the rotation is around the z axis which is per- 
pendicular to the plane of the motion. 

2.7 The gait 

A biological creature swims through a fluid by changing its shape. The oscillatory motion 
of a fin, for example, acts to propel the creature in a forward motion. Slight variations 
to the motion allow it to accelerate, decelerate, and to turn. Fish such as eels have 
a gait which is similar to a wave travelling from head to tail with moderately large 
amplitude along the entire length. A fish such as a mackerel has a gait which is similar 
to travelling wave with small amplitude until roughly half way down the body when 
it increases sharply. For the present two-dimensionsal, three-body swimmer considered 
here, the motion depends upon the angles tpi and tp2- The particular specification of these 
angles is referred to as the 'gait'. The two gaits considered here are similar to that of an 
eel. 

The forward gait of motion was specified with 

<Pi = 2 (O)-0i(O)+/3(cos(wi)-l), (2.33) 
y? 2 = 3 (O)-02(O)+/3sin(wt), (2.34) 

where for these calculations (3 = 1 and ou — 1. We choose 0i(O) = — /3, 02 (0) = 
and 03(0) = 0. Note that this specification is identical to that of Kanso et al. (2005) 
and Eldredge (2007), but the notation is different. The turning gait is the same except 
^(0) = 0, 2 (O) = and 3 (O) = -f3 and we take (3 = 1 and u = 1. With this gait the 
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angles ipj are never positive. This turning specification is the same as that of Kanso et 
al. (2005). 



2.8 The kernel 

In this paper we use the fourth order Wendland kernel for two dimensions (Wendland 
1995). With q = r/h this kernel is given by 

W(\r\,h) = ^(l + 2q)(2-q)\ (2.35) 

when q < 2 and zero otherwise. We take the initial h = 1.5dp. 

2.9 The initial conditions 

In the present simulations the liquid SPH particles were initially placed on a grid of 
squares of side dp which defines the liquid particle spacing. Those at least dp outside the 
boundary formed by the ellipses and the skin were retained, the boundary particles had 
a spacing A = dp/n where n was typically 2. The mass of the fluid particles was podp 2 , 
and the mass of the boundary particles was 1/n of this mass. 

The boundary force particles on the body were placed around each ellipse with a 
spacing as close as possible to dp/n. Having chosen which body particles connect to the 
skin, the skin sections were placed on straight lines a indicated in Figure 1. The fluid 
particles are not initially in equilibrium with the boundary forces so we allow them, and 
the skin particles, to move under damping. The rule for damping is given in the following 
section. 

After the damping is finished, the motion starts with the initial conditions set so 
that the fluid and skin particles have zero velocity and the bodies have zero net angular 
momentum and linear momentum consistent with the time derivatives of the constraints. 
The details of this are given by Kajtar and Monaghan (2008). 

2.10 The time stepping 

The time stepping is based on the second order symplectic integrator often called the 
Verlet integrator. The basic equations we integrate take the following form for a liquid 
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particle. Throughout this section, for any quantity A, A denotes its value at the begin- 
ning of the time-step, A 1 / 2 at the mid-point, and A 1 at the end of the step. The other 
particles do not change their density. 



dr a 

dt 

dt 

dp a 

dt 



= v„, (2.36) 
= F a , (2.37) 
= V a , (2.38) 



In the first stage of the integration, the mid-point values are calculated for r a , p a 
and h a , the body positions and orientations R&, 9 k , and the relative boundary particle 
positions dj = rj — R&. With 8t denoting the time step 



r l/2 


= r° a + l8tv° a 


(2.39) 






(2.40) 


hT 




(2.41) 




= R° k + lStV° k 


(2.42) 


or 


= + 


(2.43) 




= dj + ^zxdj 


(2.44) 






(2.45) 



With the mid-point coordinates known, ^J 2 , f^ 2 , F]/ 2 and r]j 2 can be calculated. The 
last two involve the Lagrange multipliers and their calculation is discussed by Kajtar and 
Monaghan (2008). 

The time-step for v and r is then completed by 

v\ = vl + StFl' 2 , (2.46) 
I = r^ + ytvl (2.47) 



With v 1 and r 1 known D\ can be calculated (this requires another sweep over the particles) 
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and the step for p and h completed according to 

p\ = P T + ¥tVl (2.48) 
hl /2 

K = TTWvUtiY (2 ' 49) 

The step for the body velocity, coordinates, angles and angular velocity is completed by 

Vi = V» + * I V fV> + F ^ | , (2.50) 



R\ = Rf + lStVl (2.51) 

Ik \jes k J 
61 = 9l /2 + l5tQl, (2.53) 

and the positions and normals of the body boundary particles at the end of the step are 
given by 



r 



3 3 



1 + 

d} + R\. (2.55) 



The damping is achieved by replacing (2.45) by 

= (v° + 5t^J 2 )D, (2.56) 



where 



The function w is given by 



D — 1 — e~ w . (2.57) 



w = Mfo-'O (2.58) 
rod 

where is the number of damping steps (typically ~ 2000), and n' is the current step. 
D is set to 1 for n' > n d . The damping steps may seem large, but for these calculations 
which involve ~ 12000 steps it is not significant. However, it would be desirable to have 
more efficient damping. 
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The time step size, St, is updated at the end of each time-step by 

x . 1 . ( h a b r ab A s \ 

St = -mm [ , — , — , 2.59 

where the minimum is over all fluid, boundary and skin particles evaluated at the mid- 
point of the time step. The first and last terms are CFL conditions for wave propagation 
in the fluid and in the skin respectively. The second term, r a b/V, ensures that St is 
sufficiently small to follow the motion of particles very close to a boundary. 

3 Numerical tests 

In the absence of skin our algorithm has been tested (Kajtar and Monaghan 2008, 2010) 
by detailed comparison against the results of experiments and those obtained by other 
authors. These include the motion of a tethered cylinder in a channel, the forced oscillation 
of a cylinder, and the inviscid calculations of Kanso et al., (2005) and Melli et al., (2006) 
where our SPH results showed convergence to the inviscid results at Reynolds numbers 
of ~ 5000. They were also in good agreement with the viscous calculations of Eldredge 
(2006, 2007, 2008). 

In this paper we begin with a test of our model of elastic skin by following the 
approach to equilibrium of fluid in a tank with an elastic skin bottom. We confirm the 
convergence of the calculations with finer resolution and show the the final displacement 
agrees with approximate theory. We then describe the simulated motion of linked ellipses 
with and without skin. 

3.1 Static tank with an elastic base 

The tank had depth D = 1 and width L — 1, while the fluid density was p = 1000kg/m 2 
(because the system is two dimensional a unit thickness in the third dimension is assumed) 
giving a total fluid mass of 1000 kg. The tension of the skin was T = 2pgD and its density 
p s = 1000kg/m 2 , and thickness i s = 0.05m. With these parameters the speed of wave 
propagation along the skin is v s ~ 19.8m/s. The Reynolds number is Re = 50. We take 
the speed of sound to be c = 10^/gD. The fluid particles were placed on a grid of squares. 
In order to determine the convergence the calculations were run for a number of different 
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initial particle spacings in the range dp — 1/20 to 1/60. The ratio of the fluid particle 
spacing to the boundary (and skin) particle spacing was 2. 

The skin was initially horizontal and when released the skin and fluid began damped 
oscillations which were followed until the skin was in equilibrium. The variation of the 
period of oscillation with the square of the initial particle spacing is shown in Figure 
3. The convergence is second order. The final position of the centre of the skin can be 
estimated by linearizing the equations of equilibrium of an elastic skin. Solving these 
equations we find that the displacement of the skin r](x) is given by 

«*) = '-4£^*-L). (3.1) 

At the highest resolution the SPH result for the skin displacement at x — L/2 is -0.068, 
which differs from the value -0.066 from (3.1) by 3 percent, which is satisfactory bearing 
in mind that (3.1) is only approximate. 

3.2 Motion of the linked bodies 

Our aim is to compare the motion of the linked bodies, with skin and without, for both 
the forward and turning gaits. In the first set of tests, the linked elliptical bodies were of 
equal size and mass. The ellipses had semi-major axis a = 0.25, semi-minor axis b = 0.2a, 
and distance between the tip of the ellipse and the pivot c = 0.2a. The second set of 
tests considered unequal sized bodies, but with the same total mass as for the first set. 
The body length parameters were a\ = 0.25, a 2 = a 3 = 0.5ai, b\ = 0.4ai, 62 = 0.5ci2, 
63 = 0.3a 3 , and for all bodies c = 0.2ai. In all cases, the densities of the bodies were the 
same as the fluid, p = 1000. 

The Reynolds number 9ft is defined by using the characteristic velocity V = 2au and 
the characteristic length L = 2a, so that 

3ft (3.2) 

In the present simulations 9ft = 200. The speed of sound was c s = 20au, and the bound- 
aries of the ellipses were defined by boundary particles with spacing as close to dp/2 as 
possible. The motion takes place in a domain with periodic rectangular cells. Based on 
the convergence studies of Kajtar and Monaghan (2010), the initial particle spacing was 
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Figure 3: The period of the damped oscillations of the skin forming the base of a tank of fluid. The 
horizontal axis is the square of the initial particle spacing dp. The convergence is clearly second order in 
dp. 
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Body 


1 


2 


3 


n 


-2.771 x KT 1 


-2.771 x lO" 1 


7.229 x lO" 1 


X 


-4.664 x 10~ 2 


2.332 x 10~ 2 


2.332 x 10~ 2 


Y 


-4.079 x 10~ 2 


8.726 x 10~ 2 


-4.647 x 10~ 2 



Tabic 1: The initial velocities for the forward gait 



Body 


1 


2 


3 


n 


-2.143 x 10" 1 


-2.143 x lO" 1 


7.857 x 10- 1 


X 


6.612 x 10~ 2 


6.612 x 10~ 2 


-1.322 x 10- 1 


Y 


-6.469 x 10~ 2 


6.388 x 10~ 2 


-8.063 x 10" 4 



Table 2: The initial velocities for the turning gait 



chosen to be dp = 1/60. The fluid spans from x min = to x max along the horizontal 
axis, and from y min = to y max in the vertical axis. For the forward gait, the domain 
was of size x max x y max = 6 x 4.5, and the initial coordinates of the centre of mass of the 
middle body were (X2,Y 2 ) = (0.4a; max , 0.6y max ). For the turning gait, the domain was 
imax x Z/max = 5x5, and the initial coordinates were (X 2 ,Y 2 ) = (0.4x max , 0.5y max ). 

The skin segments were attached to the ellipses at the mid-point between where the 
major and minor axes intercept the ellipse. In other words, the last skin particle belonging 
to a skin segment was connected to the body particle 1/8 of the distance around the ellipse 
from the major axis. Since the skin thickness was chosen as i s = 0.05, the skin segments 
had a significant mass. For example, with equal sized bodies, the total skin mass was 
approximately 50% of the total mass of the three bodies. In order to account for this 
added mass, the body masses were reduced so that the total mass of the swimmer with 
skin or without was the same. For the cases where the body masses were unequal, their 
masses were reduced such that the body densities remained equal. The initial linear 
velocities and angular momentum are given in Table 2. 

The effect of the skin on the velocity field is shown by Figures 4 and 5. The flow 
around the tips of the end points is very similar in each case, but Figure 4 shows a flow 
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1.5 2 2.5 3 3.5 

X 

Figure 4: The velocity field for the case of forward motion in the absence of skin and equal sized ellipses. 

between the gaps which influences the flow along the body. 

Table [3] shows the results for the forward gait, with skin and without, and for equal 
and unequal body sizes. Table [4] shows the results for the turning gait. For the forward 
gait, the distance travelled by the swimmer in three periods (or t = 6tt) was taken as the 
length between the initial and final positions of the centre of mass of the middle body, 
R,2. The angle of motion (in radians) was taken as the angle of this length relative to 
the x axis. For the turning gait, the amount of rotation was taken as the angle of the 
middle body relative to the x axis after three periods (or t = 6ir). Following Kajtar 
and Monaghan (2010), the average power V expended by the swimmer was computed by 
numerically integrating the following expression 

v= Ik^ ( a ^ 1)( ^ 2 " Qi) + A ^ 2)(fis " Q2) ) dt - (3 - 3) 

Figure 6 shows the zig-zag path of the centre of mass of the central body in the case 
of skin. Comparison with Figure 7 shows that the presence of skin results in a sharper 
loop at the extremes of the zig-zag, and the overall path is at a steeper angle to the 
horizontal. 
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Figure 5: The velocity field for the case of forward motion with skin and equal sized ellipses. The skin is 
shown by a dark line which can be best viewed under magnification as provided by Adobe. 



body sizes 


skin 


distance 


V 


angle 




no 


0.7038 


83.65 


-0.289 


equal 










yes 


0.7925 


74.63 


-0.299 




no 


0.3548 


24.97 


-0.525 


unequal 










yes 


0.4799 


21.54 


-0.541 



Table 3: Distance travelled by the swimmer using the forward gait, in time t — 6ir, with skin and without, 
and for equal and unequal body sizes. The power expended and angle of motion are also given. 
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Figure 6: The zig-zag path of the centre of mass of the central body in the case where there is skin and 
the ellipses are of different size. The position of the bodies is shown for the last point of the zig-zag. 



body sizes 


skin 


rotation 


V 




no 


3.7532 


142.52 


equal 








yes 


4.0487 


135.60 




no 


4.0653 


56.58 


unequal 








yes 


3.4988 


49.48 



Table 4: Amount of rotation (measured in radians) using the turning gait, in time t = 6tt. The power 
expended is also given. 
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Figure 7: The velocity field for the case of forward motion in the absence of skin. The position of the 
bodies is shown for the last point of the zig-zag. 
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For the forward gait, it is clear that the inclusion of skin improves the efficiency 
of the swimming. For equal body masses, the motion is improved by almost 13%, and 
for unequal masses by ~ 39%. We relate these improvements to the skin preventing 
fluid passing between the bodies. Although the swimmer with unequal body masses 
travels approximately half the distance in the same period of time, the power expended 
is reduced by a factor of ~ 3. In this study, the gait, [3, and u were identical. Kajtar 
and Monaghan (2010) showed that the efficiency of the swimming depends upon these 
parameters. Undoubtedly the efficiency of the swimmer with unequal body masses would 
also improve for a different set of parameters. Swimming with skin or without only 
changes the angle of motion by ~ 3%, but the mass distribution has a much larger effect. 

For the turning gait with equal body masses, the middle body rotates ~ 8% further 
when the skin is included. As with the forward gait, the power expended is slightly 
reduced. Without skin, the swimmer with unequal body masses turns ~ 8% further. 
However, with the skin included, it rotates less. 

4 Conclusions 

In this paper we have described a Lagrangian based SPH algorithm which allows us to 
simulate a system consisting of fluid, rigid body and elastic skin. The imposed change of 
angle between the bodies is included by using Lagrange multipliers. By basing the time 
stepping on a symplectic integrator the linear momentum and angular momentum of the 
system are conserved to high accuracy; for the former the errors are due to round off and 
for the latter they are dominated by the effect of the periodic boundary conditions which 
do not conserve angular momentum for a particle system. The algorithm allows us to 
treat both straight line and turning motions and bodies of arbitrary shape and relative 
size. We have applied the algorithm to three linked ellipses which may be of different 
size, and shown that the presence of the skin affects the speed, power output and turning 
capacity of the bodies. There is no difficulty in extending the code to deal with more 
bodies. 

Within the limits imposed by the algorithm being two dimensional our results sug- 
gest many interesting applications including the study of bodies swimming through a free 
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surface, or in a stratified medium, and the analysis of the hunting gaits of predators and 
the escape gaits of their prey. 



5 Appendix 

In this appendix we write out the Lagrangian for the equations of motion in the absence 
of viscous forces, and show how the variational principle with the continuity equation as 
a constraint gives the pressure forces. 

5.1 The Lagrangian 

The Lagrangian C consists of the following terms 

C = C(fluid) + C(bodies) + C(skin) + C(int), (5.1) 
where the first three terms have the form 



C( fluid) = ^2m b Q 



b 

N b 



2 v 2 b -u(p)), (5.2) 



and 



C(bodies) = 2 MkV > + 2 h ^ k) ' (5 ' 3) 

C(skin) = {\ m ° v l K ( r ™'?\ > ( 5 - 4 ) 

The fourth term involves the interactions between the particles. The first step is to note 
that the terms r j/(|r j|) can be written in terms of the gradient of a potential according 
to 

^■/(Kl) = (5.5) 



where $(|r aj -|) is a potential energy which we will denote by $ aj . The interaction part of 
the Lagrangian is then given by 

N b N s 

" V 



N b N s 

dint) = - Yl Yl m b m i®bj m b m ^ba o Kr -v ( 5 - 6 ) 



k=l jeB(k) b k=l jaeS(k) cr j 
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where in the last term the values of a and j are those for connected pairs. Following the 
usual rules the inviscid equations of motion can be worked out. Because the Lagrangian 
is invariant to translations and rotations of the coordinate system the linear and angular 
momentum are conserved. A discrete version of Kelvin's circulation theorem can also be 
deduced (Monaghan 2005). In order to work out Lagrange's equations for the rigid bodies 
it is necessary to relate the change in position of the centres of mass of the bodies and 
their angles 9 to changes in the positions of the body particles. To do this we note from 
(2.32) that for body force particle j on body k 



Because the continuity equation must be satisfied it acts as a constraint when using 
the least action principle. We consider this next. 

5.2 Least action and the Continuity equation 

In this section we consider a purely fluid dynamical problem with the Lagrangian 



b y 7 

which is to be substituted in the least action principle and varied with the continuity 
equation 



acting as a constraint. The summation over r] denotes a summation over all the particles. 
The variational principle of least action results in the equation of motion of particle a 




(5.7) 




(5.8) 




(5.9) 




(5.10) 



where S denotes a Lagrangian change. We can write 




(5.11) 



and note from the continuity equation that 




(5.12) 
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From the previous equation we deduce that 



= £(i-gV^, (5.13) 



where 5 a b is a Kronecker delta which is 1 if a = b and zero otherwise. Substitution of 
(4.13) into (4.11) then gives 

5C m a P a ^ ^m b P b 

-=- = -^^m^aW^ + ^—^VbWab. (5.14) 

0T a Pa v b Pb 

The first term is summed over all the particles whereas the second term is summed only 
over the fluid particles. However, the pressure assigned to the boundary particles and 
the skin particles is zero so we can extend the second summation over all the particles. 
Noting that V a W a b = — VbW a b we can finally write 

This gives the pressure terms on the right hand side of (2.1) where it has been split into 
separate contributions from the fluid particles, the body particles, and the skin particles. 
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